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ABSTRACT 

Measuring the relative mass contributions of luminous and dark matter in spiral galaxies is impor- 
tant for understanding their formation and evolution. The combination of a galaxy rotation curve 
and strong lensing is a powerful way to break the disk-halo degeneracy that is inherent in each of 
the methods individually. We present an analysis of the 10-image radio spiral lens B1933+503 at 
z\ = 0.755, incorporating (1) new global VLBI observations, (2) new adaptive-optics assisted K-band 
imaging, (3) new spectroscopic observations for the lens galaxy rotation curve and the source red- 
shift. We construct a three-dimensionally axisymmetric mass distribution with 3 components: an 
exponential profile for the disk, a point mass for the bulge, and an NFW profile for the halo. The 
mass model is simultaneously fitted to the kinematics and the lensing data. The NFW halo needs to 
be oblate with a flattening of a/c = 0.33^q q5 to be consistent with the radio data. This suggests that 
baryons are effective at making the halos oblate near the center. The lensing and kinematics analysis 
probe the inner ~10kpc of the galaxy, and we obtain a lower limit on the halo scale radius of 16kpc 
(95% CI). The dark matter mass fraction inside a sphere with a radius of 2.2 disk scale lengths is 
/dm, 2. 2 = 0.43^o'o9- The contribution of the disk to the total circular velocity at 2.2 disk scale lengths 
is 0.76^0 065 suggesting that the disk is marginally submaximal. The stellar mass of the disk from our 
modeling is log lo (M*/M ) = 11.06±g;?? assuming that the cold gas contributes ~20% to the total 
disk mass. In comparison to the stellar masses estimated from stellar population synthesis models, 
the stellar initial mass function of Chabrier is preferred to that of Salpeter by a probability factor of 
7.2. 

Subject headings: galaxies: halos — galaxies: individual (B1933+503) — galaxies: kinematics and 
dynamics — galaxies: spiral — gravitational lensing 



1. INTRODUCTION 

The discovery of flat rotation curves near and be- 
yond the optical edge of galaxies pr ovides strong evi- 
dence for the existence of dark matter (|Rubin et al .11 19781 
iBosmal 119781 : Ivan Albada Sancisil ll986|). Since then, 
observations of the cosmic microwave background, su- 
per novae, galaxy clusters, weak lensing, baryon acous- 
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tic oscillations, and gravitational lens time delays 
indicate that our Universe is well described by a 
model compris ed of cold dark matter and dark en- 
ergy fsee. e.g.. [Kc-matsu et al.l 1201 ll: ISuzuki et al.l 120111: 
iMantz et alljgOlOl: ISchrabback et al.l l2Q10b iBlake et al.l 
120111 : ISuvu et al.ll2Q10f ). Even though the A-CDM model 
is successful at explaining the Universe on large scales, 
the interplay between dark matter and baryons on galac- 
tic scales remains an open question. 

TV-body simulations of dark matter particles show that 
equilibrium dark matter halos have spherically averaged 
mass density profile s that are nearly universal (NFW; 
iNavarro et al.lll996|) and are typically triaxial in shape 
(e.g.. iJing fc Sutol |2QQ2[ ). The inclusion of baryons in 
simulations is challenging due to the large dynamical 
range in mass and uncertainties in the cooling and feed- 
back mechanisms. Using a sub set of the OyerW helm- 
ingly Large Simulations project ([Schave et al. 2010) that 
included va r ious p rescriptions of cooling and feedback, 
iDuffv et"aTI (|2010l ) found that the inner profile of galaxy- 
scale dark matter halos is very sensitive to the baryon 
physics. With weak stellar feedback from supernovae, 
the inner profiles tend to steepen and become more 
isothermal as a result of the high central baryon frac- 
tions that pulls the dark matter towards the center. This 
"halo contraction" is a l so found in other studies (e.g., 
iBlumenthal eT"aTl fl986t fGnedin et al.ll200i l20TTh . On 
the other hand, with strong feedback from both super- 
novae and an active galactic nucleus, the inner profile is 
very similar to the NFW profile from dark-matter only 
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simulations. Measuring the inner profiles of dark mat- 
ter halos therefore helps determine the kinds of baryonic 
processes that occur during galaxy formation and evolu- 
tion. 

Observationally, probing the inner profiles of dark mat- 
ter halos with rotation cur ve data is difficult due to the 
disk- halo degenera cy (e.g., Ivan Albada fc Sancisil 119861 : 
IDutton et al.l |2QQ5[ ). Since the rotation curve primar- 
ily depends on the total enclosed mass within spherical 
radii, a heavy disk with a light halo and a light disk 
with a heavy halo can both be fit to the rotation curve. 
The degeneracy is prominent in models where the disk 
and halos have fixed parametric forms, and is reduced in 
self-consistent models where the hal o shape changes in 
respo nse to the presence of the disk ([Amorisco fc Bertml 
120101 ). To circumvent the disk-halo degeneracy, some 
studies have assumed that the disk contributes maxi- 
mally to the circular velocity. However, studies based 
on the Tully-Fisher relation or the central vertical veloc- 
ity dispersion of disk st ars have shown that disks tend 
to be submaximal fe.^.. ICourteau fc Rix| [r999: Bo ttemal 
119931 : iBershadv et aLll2Qllh . To break the disk-halo de- 
generacy without resorting to maximal-disk assumptions, 
one needs to measure independently the relative mass 
contribution of the disk and the dark matter halo, or 
equivalently, the mass-to-light ratio (M/L) of the disk. 
Stellar population synthesis (SPS) models allows estima- 
tions of the stellar mass and hence the M/L of the disk. 
However, uncertainties in the stellar mass due to, for ex- 
ample, the unknown stellar initial mass fu nction (IMF), 
limit s the accuracy of this approach (e.g., IConrov et al.l 
l2QQ9h . Therefore, disentangling the contributions of the 
disk and the halo to the rotation curve is key in under- 
standing both the inner halo profile and the stellar IMF. 

An effective way to overcome the disk-halo degener- 
acy is to combine rotation curves with strong gravita- 
tional lensing. If a spiral galaxy lies along the line of 
sight between the observer and a background source, the 
source can be strongly lensed into multiple images by 
the spiral galaxy (e.g., iTreul l2QlQh . While kinematics 
probe mass within spheres, strong lensing probes mass 
enclosed within cylinders (within the "Einstein radius", 
which is roughly the radial distance of the multiple im- 
ages from the lens galaxy center). The combination of 
the two methods with different mass dependence breaks 
the disk-halo degeneracy. Spectroscopic and imaging 
surveys in recent years have substantially enlarged the 
sample of spiral lenses, totaling more than 20 now (e.g . , 
iFeron et al.ll200l iSvgnet et al.l 1201(1 iTreu et al.ll2Qllh . 
The first few analyses of spiral lenses are already inform- 
ing us about disk- maximality and the stellar IMFs in 
these systems (e.g., [R oopman set al.l 119981 : iMaller et al.l 
I2000t iTrott et alJl2010flDutton et al.ll201l[ ). 

In this paper, we study B1933+503, a spiral gravita- 
tional lens at z\ = 0.7 55 with 10 radio l ensed images 
that was discovered by ISvkes et"aTl (1 19981) in the Cos- 
mic Lens All-Sky Survey (CLASS; iMvers et al.l 120031 : 
IBrowne et al.l l2003t) . Pre vious modeling of the radio 
data by ICohn et all (|200lf ) tested power-law models for 
the combined dark and baryonic mass distribution of 
the lens. Here, we obtain new radio, infrared and 
spectroscopic observations, and construct a 3-component 
mass model (for the disk, bulge and dark matter halo) 
that is simultaneously fitted to both the kinematics and 



lensing data. The an alysis is very simila r in spirit to 
the one presented by IDutton et al.l (|201l[ ) on the lens 
SDSS J2141-0001, except we use an NFW instead of an 
isothermal profile to describe the dark matter halo. The 
aims of our study are (1) measure the inner shape and 
profile of the dark matter halo, (2) determine the relative 
contributions of the disk, bulge and halo, and (3) place 
constraints on the stellar IMF. 

The paper is organized as follows. In Section [2j we 
present observations of the radio lens B 1933+503. The 
alignment of the radio and near-infrared images is de- 
scribed in Section [3j and the lens light profile in the 
near-infrared image is measured in Section EJ The 3- 
component mass model is outlined in Section [5j and the 
statistical framework for the analysis is described in Sec- 
tion [6] We present the kinematics-only and lensing-only 
results in Sections [71 and [5J respectively. We discuss the 
results and implications of the joint kinematics and lens- 
ing analysis in Section [9j before concluding in Section 

[ml 

Throughout the paper, we assume a flat A-CDM cos- 
mology with H ( ) = 70 km s _1 Mpc -1 , Q\ = 1 — = 
0.73. In this cosmology, 1" corresponds to 7.5 kpc at the 
lens redshift and 8.7 kpc at the source redshift, which is 
measured in Section 12.51 Images of the lens system are 
north up and east left. Parameter constraints are given 
as the median values with uncertainties given by the 16th 
and 84th percentiles (corresponding to 68% credible in- 
tervals (CI)) unless otherwise stated. 

2. OBSERVATIONS 

We obtained both lensing and kinematic observations 
of B 1933+503 for constraining the lens mass distribu- 
tion. We present the global VLBI observations of the 
lensed radio source in Section f2TT| near-infrared imaging 
of the lens system in Sections 12.21 and 12. 3\ and spectro- 
scopic data sets for obtaining the rotation curve of the 
lens galaxy and the source redshift in Sections 12.41 and 
12.51 We describe the archival Hubble Space Telescope 
(HST) images in Section f2~6l 

2.1. Radio observations 

We observed B 1933+503 with the global Very Long 
Baseline Interferometry (VLBI) array on 1998 Novem- 
ber 27 at 1.7 GHz with a bandwidth of 16 MHz. We 
used 17 telescopes with 10 from the Very Long Baseline 
Array (VLBA) and 7 from the European VLBI Network 
(EVN). We adopt the center of component 4 of the lensed 
images as the phase center for the observations. Observa- 
tions were conducted on a cycle of 6.5 minutes, with 1.5 
minutes on a phase calibrator (B1954+51=J1955+5131) 
followed by 5 minutes on the target source B1933+503. 
The exception was the Lovell telescope, which has a 
slower slew rate and for which every other observation 
of the phase calibrator was omitted, yielding a 1.5- 
minute+11.5-minute cycle on the phase reference and 
target. The total observing time was 9 hours, providing 
the good u-v coverage required for high dynamic range 
imaging and image fidelity. A single observation of 3C345 
was obtained for fringe finding and flux calibration. 

We reduce the data with the Astronomical Image Pro- 
cessing System (ALPS 14 ) package using standard proce- 

14 developed by the National Radio Astronomy Observatory 
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Fig. 1. Global VLBI observations of B1933+503. Center: 1.7 
GHz M ERLIN observations of B1933+503 taken from lSvkes et al.l 
(1998). Contours are plotted at -3, 3, 6, 12, 24, 48 times the rms 
noise level of 150 fiJy beam -1 . The beamsize is 139 x 113 mas at 
a position angle of —13.8°. In each of the insets, the left and right 
panels show the 1.7 GHz global VLBI observations with natural 
weighting and with a 10 MA taper, respectively. In the left (right) 
panels, contours are plotted at —3, 3, 6, 12, 24, 48, 96 times the 
rms noise level of 50 (55) fiJy beam -1 , and the beamsize is 5.7 x 2.7 
mas (20.1 x 16.5 mas) at a position angle of -3.22° (60.55°). Each 
of the global VLBI panels covers a 300 mas by 300 mas area with 
the tick marks separated by 50 mas. 

dures. The images are iteratively CLEANed and self- 
calibrated (phase-only), before a single amplitude self- 
calibration solution is performed, to remove residual 
phase and amplitude errors. The final maps are pro- 
duced using natural- weighting of the visibilities and with 
a 10 MA taper to increase the sensitivity to the extended 
emission. The natural- weighting image has a root-mean- 
square (rms) noise level of 50 /iJybeam -1 with a beam 
size of 5.7 x 2.7 mas at a position angle (PA) of —3.22°. 
The tapered image has a rms of 55 /iJybeam -1 with a 
beam size of 20.1 x 16.5 mas at a PA of 60.55°. Fig- 
ure [T] shows the final maps for each of the components in 
B1933+503, except for component 6 that is not detected. 

The radio source being lensed has a compact core with 
two extended radi o lobes on opposite sides of the core 
(jSvkes et al.l 119981 : iNairl 119981 ). The radio spectral in- 
dices make the identific ation of the image s of the com- 
pact core unambiguous ([Svkes et al.|[l9 98). and the flux 
densitie s of the core comp onents showed little time vari- 
ability (Big gs et al.l 12000). Components 1, 3, 4 and 6 
correspond to the images of the compact core, compo- 
nents la and 8 are the images of one of the lobes, and 
components 2, 5 and 7 correspond to the images of the 
other radio lobe. With component 2 counting as two 
(merging) images, B1933+503 is a spectacular 10-image 
radio lens. 

We use components 1 and 4 to align the global 
VLBI data to the previous r adio observa t ions (with 
VLA, MERLI N and VLBA) in iSvkes et al.l (fl998l ) and 
iMarlow et all (jl999h . Table □ lists the image positions 
for each of the components in the previous observations 




Fig. 2. NIRC2 Kp image o f B1933+503. The radio image 
positions from ICohn et al.l (|200in are overlaid. The alignment of 
the NIRC2 and radio images are described in Section^ The lensing 
arcs are associated with radio components 1, 3, 4 and 6, which are 
images of the compact core in the source. The arcs corresponding 
to components 3 and 6 are faint due to severe dust extinction in 
the plane of the spiral galaxy. 



compiled by ICohn et al.l (|2001[ ) and in the global VLBI 
data. The merging pair of images are denoted by compo- 
nents 2a and 2b. In the global VLBI data, several com- 
ponents (1, 4, 5 and 2a) have multiple intensity peaks, 
but some of their corresponding images (3 and 7) have 
only single peaks due to, e.g., scatter-broadening dur- 
ing propagation through the disk of the lens galaxy (e.g., 
IMarlow et all 119991 : iNorburvl 120021 ) . We adopt the flux 
weighted average of the peak positions as the image posi- 
tions and estimate the uncertainties from the separation 
of the peaks. We also set the positional uncertainty of 
component 7 and 2b to that of their counter image (com- 
ponent 5). For components that are spatially extended 
(la and 8), the uncertainty is set to the geometric mean 
o f the ma jor a nd minor axis of the beam. 
ICohn et al. 
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als o com pile d flux ratios from 

ISvkes et al.1 (fl^ . lNahl (fl998h . and lBiggs et "all (pOOOh . 



We do not list or use the flux ratios because they could be 
systematically biased due to scattering an d subst ructures 
in the lens. In fact, iKochanek fc Dalai (|2004f ) showed 
that the flux ratios in B 1933+503 are anomalous due to 
substructures, and the constraints in B1933+503 lead to 
a relativel y smooth, e l liptic al macro mass model. Fur- 
thermore, et al.l (|200l[ ) found that the exclusion of 
flux ratios had little effect on the main results of their 
mass models. 

2.2. Keck AO-assisted NIRC2 image 

We observed B1933+503 with the adaptive optics (AO) 
system at the Keck II telescope on 2005 July 31. We used 
the Kp filter (at 2.2 jam) on the Near Infrared Camera 2 
(NIRC2), and took 33 dithered frames with a field of view 
of 10" x 10" that are each 6 co-adds of 30 seconds. The 
dithering allows good sky subtraction in the reduction. 
The natural tip-tilt reference star of magnitude R~16A 
located at ~21'/3 southeast of B 1933+503 and a sodium 
laser guide star were used to correct for atmospheric tur- 
bulence. 
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TABLE 1 
Radio Image Positions 







Cohn et al. (2001) 
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Notes. Column 1 lists the radio components as shown in Figure ^ Components 5 and 7 appear twice, with the peaks forming a two-im age 
system , and the extended features forming a four-image system with the merging components 2a and 2b. Columns 2-5 are from Cohn et al. 
( 2001), and columns 6-9 are from the global VLBI observations in this paper. Columns 2 and 6 are the ID numbe rs for each multipl e image 
system. Columns 3 and 7 (4 and 8) are the relative right ascensions (declinations) in the coordinate system of ICohn et al.l (|2001h where 
the origin is close to the lens center. Columns 5 and 9 are the estimated positional uncertainties. 



We use an IRAF-based algorithm to subtract the sky, 
remove bad pixels and cosmic rays, and produce weight 
maps for each frame. We use the weight maps and 
the Drizzle package (|Fruchter fc Hookll2QQ2h to coadd 
the images together. We show in Figure [2] the driz- 
zled NIRC2 image of the lens system. By comparing 
the positions of the three stars in the NIRC2 image to 
the corresponding stars in previous Hubble Space Tele- 
scope (HST) observations of the system (Section 12. 6p . 
we determine the pixel scale of the NIRC2 image to be 
9.964 mas pix" 1 . The closest field star to B1933+503 
(~2'/9 away) in the combined NIRC2 image has a full 
width at half-maximum (FWHM) of ~57mas. 

2.3. Keck NIRC image 

In addition to the high-resolution NIRC2 image, we 
observed B1933+503 o n 1996 July 31 with the N ear In- 
frared Camera (NIRC; iMatthews Soiferl fl99l ) at the 
Keck Observatory. The NIRC detector is 256 pixels on a 
side, with a pixel scale of 0'/15 pix -1 . Thus, the field of 
view of the camera is 38'/ 4 on a side, i.e., larger than the 
NIRC2 narrow camera field of view. The system was ob- 
served in both the J and K bands, but the image quality 
of the J-band data was too poor to be useful. The Ri- 
band data consist of 59 exposures, each with an exposure 
time of one minute (5 co-adds of 12 s each). The data 
were processed in a standard fashion, including steps to 
subtract the dark current, flatten the images, and sub- 
tract the sky. For each of the images, the flat-field and 
sky frames were constructed from the frames observed di- 
rectly before and after the image. The processed images 
were aligned by measuring the position of a star that ap- 
peared in each frame and then were co-added to produce 
the final image. The photometric zero point of the NIRC 
data was determined through a c omparison to the Tw o 
Micron All Sky Survey (2MASS; LSkrutskie et "all 120061 ) . 
which has K-band magnitudes for the two brightest stars 
in the image. The uncertainty in the photometric zero 
point is estimated to be 0.15 by comparing the magni- 
tudes of the two stars. This incorporates the uncertainty 
due to the difference in throughput of the 2MASS and 



NIRC K-band filters, which we estimate to be < 0.03 by 
comparing the magnitude difference between the filters 
for several stars. 

2.4. Keck ESI spectroscopy 

We observed B1933+ 503 using the Eche lle Spectro- 
graph and Imager (ESI; ISheinis et al.ll2002D instrument 
at the Keck Observatory on 2002 June 6. The slit of 
width l'/25 was aligned along the major axis of the 
lens galaxy (PA=138°). The standard echellette mode 
yielded a spectral resolution with Gaussian width of 
30kms _1 . Five exposures of 1800s each were taken un- 
der good seeing conditions of 0'/6. 

We reduce the data using the s oftware EASI2D d evel- 
oped by D. J. Sand and T. Treu (jSand et al.l 120041 ). We 
see prominent emission lines H/3 (A4861), [O n] (A3726, 
A3729), and [O ill] (A4959, A5007) from the lens galaxy 
in the spectra. For each line, we set the center of the 
galaxy to the brightest pixel in the 2-dimensional spec- 
trum since measuring the center to much less than a pixel 
is difficult due to the presence of seeing. We then bin the 
spectrum in the spatial direction by a factor of 3 to in- 
crease the signal. Thus, the center of the galaxy should 
be well within the central binned spatial pixel of size 
~ 0'/5. Any small offsets in the centroid are taken into 
account in the modeling in Section 17^31 The systemic ve- 
locity is difficult to measure from the data directly and 
is determined in the modeling (Section 17. 3p . Figure [3] is 
the rotation curve based on these lines. 

2.5. Keck NIRSPEC spectroscopy 

We carried out infrared spectroscopic observations 
of the B1933+503 lensed source on 2005 September 
9 with the Near Infrared Spectrometer (NIRSPEC; 
iMcLean et al.lfl998[ ) on the Keck II telescope. The data 
were taken through the NIRSPEC 6 and 7 filters (N6 
and N7, respectively), which gave an approximate wave- 
length coverage of 1.56 to 2.02 /xm and 2.08 to 2.52 /mi, 
respectively. Each observation consisted of 4 x 300 s 
exposures, that used dithering along the slit to improve 
the removal of the sky background during the reduction 
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Fig. 3. — Rotation curve based on the observed emissions lines 
H/3, [O n], and [O in] in the Keck ESI observations. 
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Fig. 4. — Spectra of the lensed source from NIRSPEC on Keck 
II in the N6 (left) and the N7 (right) filters are shown by the solid 
curves in the bottom panel. The rms noise spectra are the dotted 
curves. The spectra are smoothed using a 7-pixel moving average 
with each point being weighted by the inverse variance associated 
with it. The atmospheric transmission is shown in the top panel. 
In each of the N6 and N7 spectra, five emission lines (the first [N n] 
line, Ha, the second [N n] line and the two [S il] lines) are marked. 
The set in N6 is for a source redshift of z s = 1.71, and the set in 
N7 is for z s = 2.62 (with the broad line in N6 corresponding to 
H/3). The absence of emission lines in N7 rules out the previously 
identified z s = 2.62. 

stage. The total exposure time was 1 hour through each 
filter. The slit was put at a position angle of 226.7° to 
cover the two strong infrared components from the lensed 
source. The standard star HD162208 was also observed 
four times each in N6 and N7 for calibration. 

We reduce the data with a Python-based pipeline that 
removes cosmic rays, subtracts the sky, wavelength cal- 
ibrates using the atmospheric sky-lines and extracts a 
one-dimensional spectrum. Furthermore, we use the 
standard star HD162208 to correct for the response of 
the spectrograph. The reduced spectra for both filters 
are presented in the bottom panel of Figure HJ and the 
atmospheric transmission is shown in the top panel. The 
spectra are smoothed using a 7-pixel moving average with 



each point being weighted by the inverse variance associ- 
ated with it. We identify a strong emission line in the N6 
spectrum (bottom panel, left), but do not see any spec- 
tral features in the N7 spectrum (bottom panel, right). 
Therefore, we believe that the detected emission line is 
likely Ha blended with two [N n] emission lines, corre- 
sponding to a source redshift of z s = 1.71 ± 0.01. The 
uncertainty in the redshift is conservative and accounts 
for the blending of the lines and also the contamination 
by atmospheric lines on both sides of the emission line. 

We note that the strong e mission line detecte d here has 
been previously reported bv lBiggs et al.l (|200Q[ ). based on 
an unpublished spectrum that was taken with the United 
Kingdom Infrared Telescope (UKIRT). Then the line was 
interpreted to be H/3 at redshift 2.62, due to a second 
spectral feature that was believed to be Ha in the Ri- 
band. Our much better spectral resolution and higher 
sensitivity data do not detect the second emission line in 
N7, which would certainly have been detected if it were 
Ha given the relative flux of the supposed H/3 line in N6 
and the noise level. We therefore rule out a redshift of 
2.62 for the lensed source. 

2.6. Archival HST images 

Archival HST images of B1933+503 in the F160W 
filter with the Near Infrared Camera and Multi- 
Object Spectrometer (NICMOS; Proposal ID: 9744; PI: 
Kochanek) and the F555W and F814W filters with the 
Wide Field and Planetary Camera 2 (WFPC2; Proposal 
ID: 9133; PI: Falco) are available. We used MultiDriz- 
zle 15 to combine the exposures in each filter and cor- 
rect for geometric distortion. The F555W data were dis- 
carded due to the low signal-to- noise ratio (SNR). We 
oversampled the F814W and F160W images with Mul- 
tiDrizzle to the same pixel scale as the NIRC2 image, 
and show the color image composed of the three filters 
in Figure The images of the lensed source in the plane 
of the lens galaxy (corresponding to radio components 3 
and 6) suffer from dust extinction. 

We use the stars in the WFPC2 F814W images to ob- 
tain the pixel scale of the NIRC2 image in Section 12.21 
and the alignment of radio and NIRC2 images in Section 
El Furthermore, we use the F814W and F160W pho- 
tometries in Section 14.31 to estimate the stellar mass of 
the lens galaxy (that incorporates the effects of dust) in 
Section KOI 

3. RADIO AND NIRC2 IMAGE ALIGNMENT 

In order to use both the radio image positions of the 
source and the NIRC2 image of the lens galaxy to con- 
strain the lens mass distribution, we need to align the 
radio and the NIRC2 images. We assume that the two 
coordinate systems differ by a rotation and a translation. 
By aligning the stars in the NIRC2 images to the corre- 
sponding stars in the WFPC2 images with WCS infor- 
mation, we determine the north direction on the NIRC2 
image, and consequently, the rotation between the radio 
and NIRC2 images. To determine the translational offset 
between the radio and NIRC2 images, we use the centroid 
positions of the two prominent arcs in the NIRC2 image 
which we denote as P\ and P2. The separation between 

15 MultiDrizzle is a product of the Space Telescope Science 
Institute, which is operated by AURA for NASA. 
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Fig. 5. — Color image from a RGB composite of the WFPC2 
F814W, NICMOS F160W and NIRC2 Kp images. The lensed arcs 
correspond to radio image components 1, 3, 4, and 6 of the compact 
core. Component 3 is barely visible and component 6 is reddened 
due to dust extinction in the plane of the spiral galaxy. 



Pi and P2 and the angle of the segment connecting Pi 
and P2 match those of components 1 and 4 within 5 mas 
and 0.3°, respectively. The matching of Pi and P2 to 
components 1 and 4 agrees with p revious identif i cation s 
of features in the F160W image of iMarlow et al.l (|l999f ). 
We construct the coordinate translation from the NIRC2 
to the radio image by mapping the midpoint of the seg- 
ment connecting Pi and P2 to the midpoint of the seg- 
ment connecting components 1 and 4. We show in Figure 
[2] the superposition of the radio and NIRC2 data sets. 

The positions Pi and P2 of the arcs in the Keck infrared 
image can be determined with an accuracy of 5 mas. We 
add this systematic alignment uncertainty to the posi- 
tional uncertainty of the galaxy centroid in the mass 
modeling. We also incorporate the uncertainty in the 
rotation between the NIRC2 and radio image (0.3°) to 
the uncertainty in the PA of the lens galaxy. 

Having determined the transformation from the 
NIRC2 image to the radio observations, all coordinates 
values are reported with respect to the system used in 
Table [1] for the remainder of the paper. 

4. LENS GALAXY LIGHT DISTRIBUTION 

In this section we describe the steps taken to model 
the light profile of the lens galaxy in the NIRC2 image. 

4.1. PSF 

A model of the point spread function (PSF) is needed 
to extract the intrinsic light profile of the lens galaxy 
without atmospheric and instrumental blurring. We use 
the star that is ~2'/9 northwest of B 1933+503 in the ob- 
served field as a model of the PSF. Previous works have 
shown that field stars serve as good PS F models, espe- 
cially for spatially extended objects (e.g.. iMar shall et al.l 
[20071 : ISuvu eTHll2QQ9l) . 

4.2. Lens light profile 



To mo del the lens gala xy light, we use the Galfit 
package ([Peng et al. 2002). We mask out the lensing arcs 
and the spiral arm-like features in the southeast region 
of the galaxy. Optionally, we also mask out the central 
region of the galaxy. We find that the galaxy is well 
described by two exponential disk profiles with a com- 
mon centroid (or a single exponential disk profile, if the 
central region was masked out). Specifically, we employ 
Sersic profiles with index n ser sic = 1 that correspond to 
exponential disks and also fit a uniform background for 
the sky. One of the Sersic profiles corresponds to the disk 
of the galaxy, and the other is centrally localized with a 
small effective radius (~0"05) 16 . We identify this latter 
component as a bulge, although it could also be a bar 
given the limited resolution. We estimate the uncertain- 
ties on the Sersic parameters based on differences in the 
best-fit parameter values for different choices of masks. 
This systematic uncertainty dominates the statistical un- 
certainty of the fit. For the galaxy centroid, we include 
the systematic uncertainty of 5 mas (from the alignment 
of NIRC2 and radio data) which dominates the overall 
positional uncertainty. Table [2] lists the best-fit values 
for the mask containing the lensed arcs and spiral arm 
features. 

We see in Figure [6] that the two-component Sersic 
model reproduces the overall structure of the lens galaxy 
light, and the background fit yields uniform sky residu- 
als. The reduced y 2 is 1.03 in the fitting region. The 
residuals show, apart from the three strong lensing arcs 
(corresponding to radio components 1, 4, and 6 in Figure 
[2j), some small-scale features that could be spiral arms or 
tidal features in the lens galaxy. In Sections 17.31 and 18. 2\ 
we account for the effects of these residual features in the 
mass modeling by adding systematic uncertainties to the 
line-of-sight velocities and the radio image positions. 

4.3. Lens photometry 

The photometry of the lens is required for estimating 
the stellar mass from SPS models in Section [9T4l To ob- 
tain the integrated K-band magnitudes for the disk and 
the bulge based on the exponential profiles in the previ- 
ous section, we first determine the photometric zero point 
of the NIRC2 image by calibrating it with the low resolu- 
tion NIRC K-band image. We then integrate the model 
light profiles and list in Table [3] the K-band magnitudes 
of the disk and the bulge, where the estimated uncer- 
tainty comes from calibration, difference in throughput 
of NIRC2 Kp and NIRC K filters, and variations between 
different arc masks. 

The HST images have significantly lower SNR than the 
NIRC2 image. Therefore, we use the structural param- 
eters (axis ratio, position angle, and effective radius) of 
the exponential profiles from the NIRC2 image (listed in 
Table [2]) to obtain the integrated magnitudes in F814W 
and F160W with Galfit. We use the nearest field star 
to approximate the PSF. The separation of the arcs and 
the lens light is more difficult in these low SNR images, 
and we conservatively adopt an uncertainty of 0.3 mag- 
nitudes from various choices in the arc masks. The inte- 
grated magnitudes for the disk and the bulge are listed 

16 The effective radius of the central component remains small 
(< (y/07) even when the Sersic index is allowed to vary between 1 
and 4. 
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TABLE 2 
Lens galaxy light 



AR.A. 


ADec 


Re 


^sersic Q 




(arcsec) 


(arcsec) 


(arcsec) 




(°) 


Disk 0.040 ± 0.005 


-0.036 ±0.005 


0.85 ±0.05 


= 1 0.63 ±0.06 


138.0 ±1.5 


Bulge 0.040 ± 0.005 


-0.036 ±0.005 


0.055 ± 0.006 


= 1 0.41 ±0.02 


147 ±3 



Notes. Sersic profile parameters for the galaxy disk and b ulge based on the N IRC2 image. Columns 2 and 3 are the relative right ascension 
and declination, respectively, in the coordinate system of ICohn et al.I (|2001h where the origin is close to the lens center. Column 4, 5, 6 
and 7 are the effective radius, Sersic index, axis ratio and position angle of the Sersic profile, respectively. 




(c) 













Fig. 6. — Lens galaxy light, (a) NIRC2 Kp image, (b) modeled lens galaxy light based on two exponential disk profiles for the galaxy 
disk and bulge, (c) residual image. 



TABLE 3 
Lens Photometry 



WFPC2 F814W NICMOS F160W NIRC2 K 



Disk 
Bulge 



19.0 ±0.3 
23.0 ±0.3 



17.5 ±0.3 
22.1 ±0.3 



18.5 ±0.3 
22.9 ±0.3 



Notes. Magnitudes are in the AB system. 



in Table O 

5. GALAXY MASS COMPONENTS 

We decompose the spiral galaxy into three mass com- 
ponents: a disk of stars and gas, a bulge and a dark 
matter halo. We now briefly describe the mass distribu- 
tion for each component. 

5.1. Disk 

We assume that the disk of stars and gas in the galaxy 
is thin and circularly symmetric. Furthermore, we as- 
sume that there is a constant M/L throughout the disk. 
Hence, the exponential profile for the light in Section 
14.21 implies that the profile for the elliptical surface mass 
density of the projected (inclined) disk is 

E£(ii') = £ d ,o exp[-i?7# d ], 

where £d,o is the normalization of the disk (set by the 
M/L), R d is the scale radius of the disk, and R is given 
in terms of the coordinates x' , y' along the major /minor 
axes centered on the galaxy by 



R'- 



.,'2 



(2/7<Zd) 2 



The axis ratio q d is related to the inclination angle i of 
the disk by 

<? d = cosz, (3) 

where i = 0° corresponds to a face-on disk. 

In terms of the radial coordinates in the plane of the 
disk, R = \J x 2 + y 2 , the surface mass density is 

E d (R) = E d , qa exp [-R/R d ] , (4) 
where the extra factor of q d m the normalization en- 
sures that an inclined disk and a face-on disk have the 
same total mass. The scale radius R d is related to the 
effective radius (the half-light radius in Table [2j which 
is also the half-mass radius with a constant M/L) by 
R d = Re/1.678. 

For gravitational lensing, the quantity of interest is the 
dimensionless surface mass density, 

K d (R') = ^(R')/Z crit (5) 
= ft d ,o exp[-R'/Rd] , (6) 

where 

c 2 D s 

Scrit= 4^fe^fc' (7) 
^d,o is the disk strength, c\ is the speed of light and G 
is the gravitational constant. The distances D dl D s , and 
D ds are angular diameter distances to the lens, to the 
source, and between the lens and the source, respectively. 

(1) 5.2. Bulge 

The galaxy light fitting suggests that we can model 
the bulge as a point mass, since the effective radius of the 
corresponding Sersic profile is very small compared to the 
radial range of positions spanned by the lensed images 
(~(y/3 to Of! 8) and most of the points in the rotation 

(2) curve. 



8 



Suyu et al. 



5.3. Dark matter halo 

We assume a NFW profile dNavarro et a l. 1996) with 
a triaxial shape ([Jing fc Sutoll2002[ ) for the dark matter 
halo. The three-dimensional density distribution is given 

by 

Pb( r ) = 7Z7Z wi , Z7Z ^2 ( 8 ) 



(r/r h ,o)(l + r/r h , ) 2 



where 



2 2 

r — c 



b 2 



a <b < c. 



(9) 



The parameters a, 6, and c describe the triaxial shape 
of the halo. The orientation of the dark matter halo as 
seen by a distant ob server can be descr ibed by two angles 
# and (p. Following lOguri et al.l ((2003), we choose ($,(/?) 
as the polar angle of the observer's line-of-sight direction 
in the halo coordi nate system (x, y, z) (see Figure 1 of 
lOguri et all (120031)). W ith this definition, equations (1) 
and (2) of [Og uri et al.l (|2003l ) relate the coordinates in 
the frame of the halo (x, y, z) to the coordinates of the 
distant observer {x',y',z'). 

The two dimensional surface mass density of the pro- 
jected triaxial halo is elliptical, and the axis ratio (%) 
and position angle (</>h) of the ellips e can be calculate d 
given the values for a/c, b/c, $ and cp (|Oguri et al.ll2QQ3h . 
We use the Einstein radius for the corresponding spher- 
ical mass model, i?h,E, to characterize the strength of 
the halo since lensing can robustly measure this quan- 
tity. The Einstein radius is the radius of the ring that is 
formed by a point source lying perfectly behind a spher- 
ical halo with strength i?h,E- The normalization ph,o is 
related to the Einstein radius i?h,E by 



l-log(2)P if R h)E = l 

1 



Ph,o A < 



arctanh 
log(V) 
arctan 



l-Rh,E 
l+Rh,E 

2 

y 7 ^ 



if R h , E < 1 



Rh,E~ i 
Rh,E + l 



if Rh,E > 1 



where 



A - Ecrit 52 

4r h ,oAi ' 



(10) 
(11) 



and Rh,E = R^e/^o- Both i?h,E and rh,o are in arcsec- 
onds (which can be easily converted to, e.g., kpc, with 
the angular diameter distance to the lens, D^). 

5.4. External shear 

For the lensing analysis in Section [5J we also include an 
external shear component with strength 7 ex t and position 
angle ex t- A shear angle of ext = corresponds to an 
elongation of the images in the east-west direction. 

5.5. Combined mass model 

For simplicity, we impose a symmetry condition on the 
model. To be precise, we require that the centroids of the 
three components coincide and that the total mass dis- 
tribution is three-dimensionally axisymmetric. In partic- 
ular, the NFW profile is either prolate (a = b) or oblate 



TABLE 4 
Priors on model parameters 



Description 



Parameter Prior 



Centroid in 6\ 0\^ c 

Centroid in 62 #2,c 
Disk axis ratio 

Disk strength Kd,o 
Disk scale radius 

Bulge Mass M b 

Halo flattening a/c = 

Halo orientation angle 1? 

Halo orientation angle cp 

Halo Einstein radius i?h,E 

Halo scale radius rh,o 

External shear strength 7 ex t 

External shear angle e xt 



G(-0'/ 040, 0'/ 005) 
G(-0"036, 0"005) 
G(0.63,0.06) 
[7(0, 00) 
G(0'/51,0'/03) 
£7(0, 00) 
a/b £7(0.25,1) 
G(50°,5°) 
G(145°,3°) 
£7(0, 00) 

£7(0 / /l, 8 r/ ) = £7(0.75, 60) kpc 
£7(0,0.3) 

£7(0, 2tt) 



Notes. G(/i, a) is a Gaussian distribution with mean 11 and stan- 
dard devation a. U (a, b) is a uniform distribution between bound- 
aries a and b. In cases where b = 00, the upper boundary is set to 
a real number that corresponds to masses > 1O 13 M0, i.e., larger 
than galaxy-scale masses. 



(b = c). Based on the symmetry assumptions and the 
axis ratio and PA of the projected disk (in Table [2]), the 
o rientation of the ha lo is (tf, <p) = (50 ± 5°, 145 ± 3°). 

ICohn et al.l ((2001) found that the orientation of the 
quadrupole moment of their lens model (one-component 
total mass profile and external shear) agrees with the 
position angle of the lens galaxy within ~5°. This im- 
plies that the position angle of the projected total mass 
distribution is aligned with the light distribution of the 
galaxy. We confirm the alignment by modeling the lens 
syste m using a pseudoisother mal elliptic mass distribu- 
tion (Kassiola & Kovner 1993) in the presence of external 
shear. The alignment of the projected total mass distri- 
bution and the light implies that the NFW halo is oblate 
(in the inner region probed by lensing) since a prolate 
halo would lead to a ~90° difference in the position an- 
gles of the projected mass and of the light distributions. 

We impose Gaussian priors on (i) the centroid of the 
total mass distribution, (ii) the projected axis ratio and 
scale radius of the disk, and (iii) the orientation of the 
NFW halo, based on the observed light profile in Table 
[21 The position angle of the projected disk is set by the 
orientation of the NFW halo (based on the axisymmetry 
assumption) . We impose uniform priors on the remaining 
parameters, and summarize the priors in Table HI While 
all these priors are imposed for the lens modeling, some 
are not needed for the kinematics, such as the external 
shear 17 . In total, we have 13 mass parameters: 6 with 
Gaussian priors, and 7 "free" parameters with uniform 
priors. 

6. BAYESIAN ANALYSIS 

We use Bayesian analysis to infer the mass model pa- 
rameters. In particular, we sample the posterior proba- 
bility distribution function (PDF) of the 13 mass param- 

17 B1933+503 is not in any obvious galaxy group or cluster, so 
any external shear is likely to come from structures that are not 
dynamically associated with the lens system. 
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eters rj which is given by 



P(rj\d) 



likelihood prior 

p{d) 



(12) 



evidence 

where d is the data. The expressions for the likelihoods 
of the kinematics and the lensing data are in Sections 17.31 
and !8.2| respectively. The prior P(rj) is given by 



13 



P(r 1 ) = l[P(r ll ), 



(13) 



where rji is the i-th parameter, and P(r]i) is either a Gaus- 
sian or a uniform distribution as described in Section [531 
The Bayesian evidence is used for comparing different 
forms of models, which does not concern us in this paper 
since we consider only one form of mass model: an ex- 
ponential disk, a point mass bulge and an oblate NFW 
halo. 

7. ROTATION CURVE MODELING 

In this section we describe the modeling of the galaxy 
mass distribution based on the rotation curve con- 
straints. We assume that the gas in the disk (from which 
we observed the emission lines for the rotation curve) is 
in circular orbits, and use the circular velocities of the 
mass model to predict the line-of-sight velocities for con- 
structing the likelihood of the kinematics data. We then 
sample the posterior PDF of the kinematics data. 

7.1. Circular velocities 

For a test mass in the plane of the disk at a radius 
R from the center, its circular velocity v c (R) has three 
contributions: 

^ = ^b + ^d + ^ (14) 
where v D ,^d, and are the circular velocities of the 
bulge, disk and halo mass distributions, respectively. 

For the point-mass bulge, the rotational velocity con- 
tribution v D takes on a simple form 

GMh d5) 



R 



where G is the gravitational constant and Mb is the mass 
of the bulge. 

The contribution of the disk also has an analytic ex- 
press ion in terms of Bessel func tions Jo, ^o, A, and K\ 
(e.g., iBinnev fe Tremaindll987[ ): 

v\ = 4wGi: c R d y 2 (I (y)K Q (y) - h^K^y)), (16) 

where E c is the central surface mass density of the disk 
(= ftd,o^crit#d) 7 ^d is the scale radius of the disk, y is 
related to the radial distance R by y = t^ - , and R e is 

the effective radius of the mass prof ile of the disk. 

For the halo component, we follow iBinney fc Tremainel 
(|l987f ) and integrate the oblate spheroidal three- 
dimensional mass density p(m) with m 2 = x 2 /q 2 l +y 2 + z 2 
to obtain the circular velocity: 



AirGq 



Jo 



p(m) 



zdra, 



(17) 



- ™ 2 (i - <£) 

where = a/c is the flattening of the dark matter halo. 
For the NFW mass density (equation (|8))), we numeri- 
cally compute the integral in equation (fT7|) . 



7.2. Predicted rotation curve 

To model the rotatio n cur ve, we follow 
Ivan der Marel fc van Dokkuml ([20 07) to obtain the 
predicted line-of-sight velocities from the circular veloci- 
ties which incorporate the effects of the disk inclination, 
slit and seeing. 

The line-of-sight velocity in the absence of the slit and 
seeing at a point (x\ y') on the sky (with (x f , y') = (0, 0) 
at the center of the galaxy and x' along the major axis 
of the lens galaxy) is 



vios(z' \y f ) = ~b sm(i)v c (R), 



(18) 



where R = (x, y) is the corresponding coordinate in the 

plane of the disk with x = x', y = and R = \R\. 

We weight v\ os by the modeled galaxy light, and con- 
volve with both a Gaussian PSF with FWHM of 0'/6 and 
a square top-hat function of size W x P (for slit of width 
W and pixel scale P) for each (binned) spatial pixel along 
the slit. The resulting weighted and convolved velocity 
is the predicted line-of-sight velocity ^ r g ed at the corre- 
sponding spatial pixel. 

7.3. Kinematics analysis and results 

The uncertainties of the data points in the rotation 
curve of Figure [3] are only statistical. We include a sys- 
tematic uncertainty of 40kms _1 (which is comparable 
to the statistical uncertainty) to account for deviations 
from our assumptions of the smooth axisymmetric mass 
model, slight offset (if any) of the galaxy centroid in the 
rotation curve, and the unknown systemic velocity. We 
explore a range of values for the systemic velocity, and 
choose the value that optimizes the rotation curve fitting. 
We add the statistical and the systematic uncertainty in 
quadrature to obtain the final uncertainty on the line- 
of-sight velocity. The amount of systematic uncertainty 
is set such that the rotation curve can be modeled with 
a reduced x 2 ^l- The inclusion of the systematic uncer- 
tainty is crucial for not underestimating the uncertainty 
on the resulting mass model parameters. 

The likelihood for the rotation curve data is 



P(d B \ri) = — exp 



4 s 



,,obs 
J los,i 



pred\ 1 
los,i ) 



, (19) 



n obs 



IS 



where TVd is the number of data points (= 15), v[ os • 
the observed line-of-sight velocity of data point z, 
is the predicted line-of-sight velocity from our model and 
observational setup, oi is the uncertainty in the velocity, 
and Zd is the normalization given by 



(2^/ 2 n< 



(20) 



We use Mu ltiNest dFeroz fc Hobsonl [2008; 



iFeroz et all [2009) to sample the posterior PDF of 
the 13 model parameters. Figure [71 shows the resulting 
constraints on four of the parameters: halo flattening 
a/c, halo Einstein radius i?h,E 5 halo scale radius rh,o, and 
disk strength ftd,o- The disk-halo degeneracy appears 
as the anti-correlation between the disk strength (ftd,o) 
and the halo Einstein radius (i?h,E)' the more massive 
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halo a/c halo R h E [arcsec] halo r h [arcsec] 
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0.4 0.6 0.8 
halo a/c 



Fig. 7. — Marginalized posterior PDF for the halo flattening 
a/c, halo Einstein radius i?h,E> halo scale radius r*h o? an d disk 
strength based on only the rotation curve data. The three 

shaded areas show the 68.3%, 95.4% and 99.7% credible regions. 
The disk-halo degeneracy is illustrated in the top-middle panel: 
higher i?h,E corresponds to lower ^d,o- 

the disk, the less massive the halo, and vice versa. The 
halo scale radius is poorly constrained. Furthermore, 
the kinematics data provide little information on the 
halo shape. 

Figure [8] shows the predicted rotation curve of the 
most probable mass distribution with a reduced x 2 = 1-0 
based on the kinematics data. The bulge contributes very 
little to the mass of the system, consistent with the small 
effective radius observed for the bulge light distribution. 
At small (large) radii, the disk (dark matter halo) dom- 
inates in the enclosed mass. 

8. LENS MODELING 

In this section we discuss the properties of our three- 
component mass model based on the lensing constraints. 

8.1. Lensing deflection angles 

We briefly describe how to obtain the deflection an- 
gles for the three mass components in our model. The 
deflection angles are used to solve for the predicted ra- 
dio image positions that are needed for constructing the 
likelihood of the lensing data. 

8.1.1. Disk 

The lensing deflection angles for an elliptical mass com- 
ponent following an exponential disk profile has no ana- 
lytical form. Hence, to model the stellar disk component, 
we employ a "chameleon profile" mimicking the mass 
profile of an exponential disk whose deflection angles are 
analytical. This profile is i nspire d by the chameleon pro- 
file used in iMaller et all (|2000h . We use the following 
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Fig. 8. — Rotation curve of the most probable mass model 
based on kinematics data only. The contributions from the disk 
(blue dashed), halo (red dotted) and bulge (magenta dot-dashed) 
are indicated. Most of the mass is in the disk and the dark matter 
halo. 



approximation for the dimensionless surface mass den- 
sity of the disk: 



Kd(R') = ^d,o exp [-R' / Rd] 



*«d,0Sl 



Vr^ 



Vr^ 



(21) 
(22) 



where s 1 = 4.6849# d , s 2 = 1.1720# d and s 3 = 1.4518i? d . 
The projected enclosed mass within radius R' for the 
chameleon profile agrees with that of the exponential 
profile within 2% for the range of radii spanned by the 
lensed images (R r = 0.6Rd to 1.6i?d)« Each term in the 
square brackets in equation ([22]) is in the form of a pseu- 
doisothermal elliptic mass distributi on whose deflection 
angles can be computed analytically ([Kassiola fc Kovnerl 

Eaal. 

8.1.2. Bulge 

The bulge is modeled as a point mass , and the deflec- 
tion a ngle has a simple closed form (e.g. JSchneider et al.l 
l2006h 

«(0 = — ~T~ ( 23 ) 
c l s 

where G is the gravitational constant, c\ is the speed of 
light, Mb is the mass of the bulge, and £ is the impact 
parameter (i.e., the distance between the light ray and 
the point mass on the lens plane). 

8.1.3. Dark matter halo 

We refer the reader to lOguri et all (|2QQ3f ) for the lens- 
ing properties of the triaxial NFW halo in equations (j8j) 
and ([9]). The projected surface mass density of the halo 
is elliptical, and the deflection angles can be computed 
numerically. 

8.2. Lensing analysis and results 
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We use the radio image positions listed in Table Q] to 
constrain the mass model parameters. We add a po- 
sitional uncertainty of 10 mas in quadrature to the un- 
certainties in the table to account for systematic effects 
such as the residual features in the lens gal axy light fit 
(Sect ion 14. 2|) . presence of substructure (e.g.. IChen et al.l 
l2QQ7h and scatter-broadening through the disk of the lens 
(|Marlow et al .111 9991 : [Biggs et al.ll2004 e.g.,). The lensing 
arcs in the NIRC2 images could in principle be used in 
addition to the radio image positions to constrain the lens 
mass distribution; however, these arcs (especially the one 
associated with components 3 and 6) in practice suffer 
from dust extinction. Only the F160W image of the HST 
dat a has sufficient SN R to be used for dust correction, 
and lSuvu et a l. (2009) showed that dust correction based 
on two bands are prone to systematic effects. Therefore, 
we do not include the dust-affected NIRC2 arcs for con- 
straining the mass distribution. 

The likelihood of the radio image positions is 



P(d L \r}) = 7^ ex P 



•jobs 



zipred 
ij 



(24) 

where N sys is the number of multiply imaged systems, 
N? m is the number of multiple images in system j, 0^ 
is the observed image position, 0^ d (r]) is the modeled 
image position, cr^j is the uncertainty in the observed 
image position, and is the normalization given by 



Z L = (2tt) 



nn 

.7 = 1 i=l 



with 



(25) 



(26) 



The source position for each system of multiple images 
is needed to predict the image positions. We model the 
source position as the weighted average of the mapped 
source positions from the observed image positions and 
the deflection angles from the mass distribution. Specif- 
ically, for each image system, we take the average of the 
mapped source position j3i weighted by J~jTi/(Ti, where \i{ 
and Gi are, respectively, the modeled magnification and 
the positional uncertainty of image i. 18 In doing so, we 
have in effect marginalized the source position parame- 
ters by approximating the lensing likelihood as having a 
delta function at the weighted source position for each 
image system. 

We model the 3-component mass distribution of the 
spiral galaxy based on the lensing data with Glee 
( "Gravitational Lens Efficient Explorer" ) , a software de- 
veloped by A. Halkola & S. H. Suyu 19 . For a set of mass 
model parameter values, Glee computes the values of 

18 Our tests using analytic mass profiles show that the weighted 
source positions and the optimized source positions give consistent 
constraints on the mass parameters. Since the computation time of 
the weighted source positions is an order of magnitude shorter than 
that of the optimized source positions and the deflection angles of 
the oblate NFW halo are computationally expensive (due to the 
numerical integration), we use the weighted source positions. 

19 The software models the mass distribution in gravitational 
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Fig. 9. — Marginalized posterior PDF for the halo flattening 
a/c, halo Einstein radius i?h,E> halo scale radius rh o? an( l disk 
strength k^^o based on only the lensing data. The three shaded 
areas show the 68.3%, 95.4% and 99.7% credible regions. The 
panels are plotted on the same scales as in FigureOfor comparison. 



the lensing likelihood in equation (J24J) and the prior in 
equation (fl3|) ; these are the two ingredients needed to ob- 
tain the posterior PDF in equation ([T2]). As in the kine- 
matics analysis, we us e MultiNest (|Feroz fc Hobsonl 
l2Q08b iFeroz et al.ll2009h to sample the posterior PDF. 

Figure [9] shows the resulting constraints on the same 
parameters as those in Figure [71 The degeneracy be- 
tween the disk strength and the Einstein radius of the 
dark matter halo is visible, though not as strong as in 
the case of the kinematics-only analysis. The halo is 
highly flattened (a/c ~ 0.3). The flattening is degener- 
ate with the disk strength as shown in the top-left panel: 
a flattened halo is less massive and requires a more mas- 
sive disk. The flattening is also degenerate with the halo 
Einstein radius: massive halos with high Rh,E need to be 
more flattened to reproduce the overall ellipticity of the 
projected mass as constrained by the lensing data. The 
scale radius of the dark-matter halo is not constrained, 
as expected since lensing only probes the distribution in 
the radial range spanned by the images, i.e., ~0'/5, or 
~4kpc. Nonetheless, a small value of < 1" is rejected by 
the data at 95% CI. 

The most probable lensing model (with highest pos- 
terior PDF) has a reduced x 2 = 0.9. We show in Fig- 
ure [10] the critical curves (solid) and caustics (dashed) of 
the most probable lensing model. The open symbols are 
the observed image positions, and the solid symbols are 
the modeled source positions. The figure illustrates the 
configuration of the 10-image system in relation to its 



lenses using image positions (Halkola et al. 2008) or extended im- 
ag es with pixelated sourc e JSuvu et al. 2006). We refer the reader 
tofSuvu & Halkola ( 2010) for an example of each approach. 
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Fig. 10. — The critical curves (solid) and caustics (dashed) of 
the most probable model based on lensing data only. The open 
symbols are the observed image positions that are lab eled accord- 
mg to Table □ 

squares are the images positions in Co hn et alJ 
( 2001), and circles are the global VLBI image positions. The cor- 
responding solid symbols are the modeled source positions. The 
three groups of source positions, which have one group outside the 
astroid caustics (with image multiplicity of 2) and two groups in- 
side the caustics (with image multiplicity of 4), form the 10-image 
system. The model reproduces the observed image positions to 
within / /02, except for the merging components 2a and 2b (in sys- 
tems IV and VII) which are reproduced to within (//06 due to the 
higher positional uncertainty from the high magnification. 

3-component source (the source positions are clustered 
into three groups). The first group of sources is outside 
the astroid caustics and produces components la and 8. 
The second group is inside the astroid caustics and pro- 
duces components 1, 3, 4 and 6. The third group is near 
a fold of the caustics and produces components 2, 5 and 
7. 

9. COMBINING KINEMATICS AND LENSING 

In this section, we present results on the mass distri- 
bution for the spiral galaxy based on the kinematics and 
lensing data sets. Since the two data sets are indepen- 
dent, the likelihood is 

P(d D ,d L |r7) = P(d D \r 1 )P(d L \r ] ), (27) 

where the kinematics and lensing likelihoods on the right- 
hand side are given by equations (fT9|) and (|24|) . respec- 
tively. Figure [IT] shows the result of MultiNest sam- 
pling of the posterior. The reduced x 2 of the joint data 
set is 0.8. The marginalized parameters with 68% CI are 
listed in Table El 

9.1. Breaking the disk-halo degeneracy 

Although the kinematics constraints are significantly 
weaker than the lensing constraints, the ftd,o~^h,E panel 
(top-middle) in Figure [11] shows that the kinematics and 
lensing contours are tilted at different angles. Thus, the 
combination of the two in principle breaks the disk-halo 
degeneracy. In the case of B1933+503, most of the con- 
straints on the mass distribution come from the lensing 
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Fig. 11. — Marginalized posterior PDF for the halo flattening 
a/c, halo Einstein radius i?h,E> halo scale radius r^cb an d disk 
strength tt d based on kinematics only (red), lensing only (blue), 
and joint lensing and kinematics (black). The three shaded areas 
show the 68.3%, 95.4% and 99.7% credible regions. 

data: the 10 radio images, which both span a large range 
of radii and have positional uncertainties of only a few 
mas, provide a strong leverage in discerning contribu- 
tions from the disk, bulge and halo. This is in contrast 
to typical lenses with only two or four imag es which re- 
sult i n stronger disk-halo degeneracies (e.g.. iTrott et al.l 
l2Q10h . On the other hand, the SNR of the rotation curve 
data for B 1933+503 is only modest given the high lens 
redshift of z\ = 0.755. Nonetheless, the kinematics data 
in B 1933+503 are informative in excluding the high disk 
masses allowed by lensing, which has interesting conse- 
quences that are discussed in Section 19.41 

A disk is considered to be maximal if the fractional 
contribution to the circular velocity of th e disk at 2.2 R <\ 
is v d (2.2R d )/v tot (2.2R d ) = 0.85 + 0.1 (jSackettJ fl997h . 
Based on our lensing and kinematics model, B1933+503 
has v d (2.2R d ) = 248^6 kms_1 and v tot (2.2R d ) = 
326 ± 8kms _1 . In comparison to the Milky Way's 
circular velocity at the posit ion of the sun, y c (Ro) = 
(219±20)iV(8kpc)kms- 1 (|Reid et al.lll999f ). the spi- 
ral galaxy in B 1933+503 is significantly more mas- 
sive. The resulting disk contribution in B1933+503 of 
v d (2.2R d )/v tot (2.2R d ) = 0.76±He suggests that the 
disk is marginally submaximal. 

9.2. Shape and profile of the dark matter halo 

In the first column of Figure [TTJ we see that the dark 
matter halo in B1933+503 is oblate with a/c = 0.33t£;^ 
in order to fit to the lensing observations. The axis ra- 
tio of the projected surface mass density of the halo is 
0.67 ± 0.04, which is slightly rounder than the axis ratio 
of projected disk of 0.53 + 0.03. It appears that only pro- 
jected mass distributions of the lens with axis ratios of 
^0.6 are consistent with the lensing data. If we consider 
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TABLE 5 
Marginalized parameters (68% CI) 



disk q& disk ^d,o 


bulge M b 
(1O 1O M ) 


halo a/c 


halo gh 


halo <£h halo i?h,E 
(°) (arcsec) 


halo rh,o 
(arcsec) 


7ext 4>ext 
(°) 


0.53 ±0.03 1.1 ±0.2 


1.6 ±0.3 


33 +u - uv 


0.67 ±0.04 


137 ±2 0.4 ±0.2 


5±2 


0.02 ±0.01 103 ±14 
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Notes. The marginalized model parameters for the 3-component mass distribution constrained by lensing and kinematics. Columns 1-2 
are the projected disk axis ratio and strength. Column 3 is the bulge mass. Columns 4-8 are the halo flattening, projected axis ratio, 
projected position angle, Einstein radius and scale radius. Columns 9-10 are the external shear strength and position angle. 

the subset of the MultiNest sample with a/c ~ 0.5 
(which corresponds to a projected axis ratio of ~0.75), 
then the lensing reduced x 2 increases from 0.8 to 2.4 due 
to misfits in systems II, V, VI and VII in Table [U The 
global VLBI data with better positional accuracies and, 
in particular, the four- image systems II and VI enforce 
the high ellipticity in the projected surface mass density 
and hence the highly oblate halo. The high ellipticity is 
robust against assumptions on the dark matter halo pro- 
file. In particular, when we replace the NFW halo with 
either a singular or cored isothermal profile, the radio 
data still require a high ellipticity for the halo. 

TV-body simulations of dark matter halos indicate that 
the halos are typically triaxial and the axis ratio between 
the short and long axis (i.e., a/c in our notation) is ~0.3 
to 0.9 (95%CL) for halo masses of 1 11 - 10 12 M (e.g., 
iBett et all 120071 : IMaccio et afl l2QQ8h . Simulations with 
baryons find that the central galaxy tends to make the 
triaxial halo essen tially oblate and less flattened (e.g., 
lAbadi et al.ll2010f ). Our study of B1933+503 with the 
resulting a/c ~ 0.3 for the oblate NFW halo suggests 
that baryons are effective at making the halo oblate in 
the inner regions of the galaxy. 

Since the kinematics and lensing data probe the in- 
ner = 7.5 kpc of the galaxy, it is not surprising that 
the scale radius rh,o of the dark matter halo in our model 
(where the density transitions from p ~ r _1 to r -3 ) is not 
constrained but has a lower limit of 2'/l = 16 kpc (95% 
CI). The virial radius is difficult to extrapolate from the 
model given the uncertain scale radius. For galaxies with 
total mass of -10 12 M (applicable to B1933+503), typ- 
ical virial radii are ~250kpc. This implies a concentra- 
tion of < 16 that is consistent with concentration- mass 
relations for galax y-scale halos (e.g.. IMaccio et aTl 120081 : 
iDuffv ct "all 12010). To better constrain r^o? a rotation 
curve that extends to larger radii is needed. 

9.3. Dark matter mass fraction 

From our 3-component model of the galaxy, we can de- 
termine the fraction of dark matter as a function of radius 
by integrating the mass enclosed within spherical radii. 
Since we employ a parametrized model and the parame- 
ters of the model are constrained by the lensing and kine- 
matics data, the enclosed mass can be computed for all 
radii. At small and large radii where there is no (or low 
SNR) data, the enclosed mass is effectively extrapolated 
based on the form of the parametrized model and the 
measured values of the mass parameters. The top panel 
in Figure [12] shows the mass enclosed for the disk and 
the dark matter halo. The bulge, which is unresolved and 
modeled as a point, has mass Mb = (1.6 ±0.3) x 1O 1O M 
for all radii. In the inner 10 kpc of the galaxy, we see that 
the disk dominates in mass, and beyond that, the dark 
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/-[kpc] 

Fig. 12. — Masses and dark matter mass fractions from the 
lensing and kinematics modeling. Top panel: spherical enclosed 
mass of the disk and halo components. The bulge is modeled as a 
point mass of (1.6±0.3) x 10 10 Mq. Bottom panel: the dark matter 
mass fraction within a sphere. The dotted (dashed) vertical line is 
the scale (effective) radius of the disk. 

matter halo dominates. 

The bottom panel shows the dark matter mass fraction 
within a sphere of radius r. The rise in the dark mat- 
ter fraction is similar to the analys i s of th e spiral lens 
SDSSJ2141-0001 bv iDutton et~all (|20lTI) . The dark 
matter mass fraction within 2.2 disk scale lengths for 
B1933+503 is / DM) 2.2 = 0.43±£j$. Within the effective 

radius, the mass fraction is /DM,e = 0.37^o'.o8> which is 
consistent with the ranges of values found in previous 
lensi ng and kinem atics analyses of early-typ e galaxies 
(e.g., [Auger et al.ll2Q10at iBarnabe et al.ll201ll ). 

9.4. Implications for the IMF 

The shapes of the low- and high-mass end of the stellar 
IMF are difficult to determine observationally since low 
mass stars are intrinsically faint and high-mass stars are 
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low in abundance. In this section, we describe the esti- 
mation of t he stellar mass in B1 933+503 based on SPS 
models for IChabrierl (|2003l ) and ISalpeterl (jl955h IMFs. 
We compare these masses to the independently measured 
disk mass from lensing and kinematics to learn about the 
IMF. 

Usin g the broadband photometries in Table [3j we 
follow I Auger et al.l ([2009) to infer the stellar mass for 
the disk based on composite stellar popu lation synthe- 
sis models from lBruzual fc Charlotl ([2003 ). These mod- 
els have been e mployed in many s t udies i ncluding SPSS 
galax ies (e.g., iKauffmann et al.l 120031 : iGallazzi et al.l 
2005|) and spiral lens analyses similar to ours (e.g., 
Dutton et al. | l201l| ). Dust extinction (intrinsic) is taken 



into account in these models, and our multiband pho- 
tometries allow a good handle on dust with only a sma ll 
broadening of uncertainties (e.g.. [Bell fe de Jongll200i| ). 
Our inference of the stellar mass provides uncertain- 
ties that incorporate parameter degeneracies in the stel- 
lar population models (e.g., between age and metallic- 
ity ) . In Figure [13j we plot the inferred stellar mass in 
red-dashed (blue-dotted) assuming Chabrier (Salpeter 20 ) 
IMF. Th ese two IMFs coyer the range applicable to spiral 
galaxies (jBell fc de Jondl200lf ). 

The lensing and kinematics analysis provides an in- 
dependent measurement of the total disk mass of 
log lo (M disk /M ) = 11.1 7to. io tnat includes stars and 
gas (thin solid curve in Figure [T3|) . To extract the 
stellar mass contribution, we assume that the cold 
gas accounts for 20 ± 10% o f the total mass (e.g., 
iDutton fc van den Boschl l2QQ9h and is distributed like 
the stars. This provides an upper limit in the mass con- 
tribution of the gas to the disk (in comparison to scenar- 
ios where the gas is more extended than the stars). For 
each sample in the posterior PDF of the disk mass, we 
draw a random number / from a Gaussian distribution 
centered on 0.2 with a standard deviation of 0.1 that is 
truncated between [0,1], and subtract the fraction / from 
the total disk mass. The gas-subtracted disk mass from 
the lensing and kinematics model is log lo (M*/M ) = 

11.06™ 

and the distribution is shown by the thick 
solid curve in Figure HH 

Comparing the thick solid curve to the dashed and dot- 
ted curves in the figure, we see that our mass model of 
B1933+503 favors a Chabrier-like IMF. We can quanti- 
tatively rank the two IMFs by computing the Bayesian 
evidence which is the integral under the product of the 
PDF from lensing and kinematics (thick solid) and the 
PDF from the SPS (dashed or dotted). The ratio of 
the Bayesian evidence for Chabrier to Salpeter is 7.2; in 
other words, the probability of the IMF being Chabrier 
is 7.2 times higher than the probability of the IMF being 
Salpeter for B1933+503. 

Using the Chabrier IMF, we obtain the rest-frame 
stellar mass-to-light ratio of the disk to be M*/Ly = 
(0.6 ± O.3)M /Lv,0. By passively evolving to z = 0, 
we obtain the present day stellar mass-to-light ratio of 
M*/Ly = (1.71qJ)M0/Lv,o, in agreement with typ- 



ical^ values found for lo cal g al axies (e.g., ITrott et al.l 
ie Yen et all [20101: ICourteau fc^Rixl 119991 : 
van der Kruit fc Freeman! 120 111 and references therein). 
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We use 0.1 Mq as the lower mass limit for the Salpeter IMF. 



Fig. 13. — Comparison of the disk mass from the lensing and 
kinematics analysis, and the stellar mass from photometry and 
stellar population synthesis with Chabrier or Salpeter IMFs. A 
Chabrier-like IMF is preferred to a Salpeter-like IMF by a factor 
of 7.2 when a 20% d= 10% contribution in mass from the cold gas is 
assumed. 

Our finding of a preference toward a Chabrier-like IMF 
for the spiral galax y is i n agreement with th e results in 
IDutton et all ([20 111 ) and lFerreras et al.l (j2010f ). Nonethe- 
less, studies of m assive ellipt i cal ga l axies favor Salpeter- 
like IMFs (e.g., iTreu et alJ l2010t lAuger et all l2010ri 
Ivan Dokkum fc Conro v 2010; Spiniello et al.ll2Qllh . This 
supports a non-universal IMF for galaxies that is depen- 
dent on mass and/or Hubble type. 

10. CONCLUSIONS 

We disentangled the distributions of baryons and dark 
matter in the spiral galaxy B1933+503 by using the 
newly acquired global VLBI observations, AO-assisted 
NIRC2 imaging, rotation curve, and source redshift. We 
constructed an axisymmetric 3-component mass model 
for the galaxy with an exponential disk, a point mass 
to approximate the unresolved bulge, and a NFW dark 
matter halo. Parameters of this model were constrained 
by a combined lensing and kinematics analysis. Based 
on this study, we conclude the following 

• The image positions of the radio source span a 
range of radii and provide strong constraints on 
the 3-component mass distribution of the lens. 

• The fractional contribution of the disk to the total 
circular velocity at 2.2 is 0.76io'.o6> suggesting 
that the disk is marginally submaximal. 

• The oblate dark matter halo needs to be highly 
flattened with a/c ~ 0.3 in order to fit to the radio 
observations. 

• The lensing and kinematics data sets probe the in- 
ner ~10kpc region of the mass distribution and 
place a lower limit of 16 kpc (95% CI) for the scale 
radius of the NFW halo. 
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• The dark matter mass fraction inside a sphere in- 
creases as a function of radius. The mass fraction 
within 2.2 R d is / DM ,2.2 = 0.43±g;^ 5 and within the 
effective radius is fvM,e — 0.371q'.o8- 

• The total stellar mass of the disk based on the lens- 
ing and kinematics data sets is log lo (M*/M0) = 
11.06^0.11? assuming that the cold gas is dis- 
tributed like the stars. 

• Based on the lensing and kinematics measurement 
of the disk mass, the Chabrier IMF is preferred to 
the Salpeter IMF by a probability factor of 7.2. 

The sample of spiral lenses has been growing rapidly 
in recent yea r s tha n ks to dedica t ed su r veys (e.g., 
Marshall et all 120091 : iSvgnet et aD 120101 : iTreu et all 
20111 ). While most of the lenses do not have sources 
that are radio loud, spatially extended lensed images 
can also be used to constrain the mass distribution (e.g., 
iDutton et al.ll20lTh . The combined lensing and kinemat- 
ics modeling methods we have developed are general, and 
important insights into the interactions between baryons 
and dark matter in the formation and evolution of spi- 
ral galaxies can be derived from a sample of lenses using 
these techniques. 
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